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A UNIFORMLY ACCURATE (UA) MULTISCALE TIME INTEGRATOR 
PSEUDOSPECTRAL METHOD FOR THE DIRAC EQUATION IN THE 
NONRELATIVISTIC LIMIT REGIME* * * § 

WEIZHU BAOt, YONGYONG GAI*, XIAOWEI JIA§, AND QINGLIN TANG’ 


Abstract. We propose and rigourously analyze a multiscale time integrator Fourier pseudospectral (MTI-FP) 
method for the (linear) Dirac equation with a dimensionless parameter e € (0,1] which is inversely proportional to the 
speed of light. In the nonrelativistic limit regime, i.e. 0 < e <C 1, the solution exhibits highly oscillatory propagating 
waves with wavelength 0(£^) and 0(1) in time and space, respectively. Due to the rapid temporal oscillation, it is quite 
challenging in designing and analyzing numerical methods with uniform error bounds in £ G (0,1]. We present the MTI- 
FP method based on properly adopting a multiscale decomposition of the solution of the Dirac equation and applying 
the exponential wave integrator with appropriate numerical quadratures. By a careful study of the error propagation 
and using the energy method, we establish two independent error estimates via two different mathematical approaches 
as and where h is the mesh size, r is the time step and mo depends on the regularity of the 

solution. These two error bounds immediately imply that the MTI-FP method converges uniformly and optimally in 
space with exponential convergence rate if the solution is smooth, and uniformly in time with linear convergence rate at 
0(t) for all £ G (0,1] and optimally with quadratic convergence rate at O(t^) in the regimes when either £ = 0(1) or 
0 < £ < T. Numerical results are reported to demonstrate that our error estimates are optimal and sharp. Finally, the 
MTI-FP method is applied to study numerically the convergence rates of the solution of the Dirac equation to those of 
its limiting models when £ —>• O’*". 

Key words. Dirac equation, nonrelativistic limit regime, uniformly accurate, multiscale time integrator, exponential 
wave integrator, spectral method, error bound 


1. Introduction. Quantum mechanics and relativity theory were the two major physics discov¬ 
eries in the last century. The first successful attempt to consistently integrate these two fundamental 
theories was made by the British physicist Paul Dirac in 1928 nsun], resulting in the equation known 
as the Dirac equation. It describes the evolution of spin-1/2 massive particles, such as electrons and 
quarks. It is a relativistic version of the Schrodinger equation for quantum mechanics, consistent with 
Albert Einstein’s special relativity. Dirac’s theory led to the rigorous explanation for the fine structure 
of the hydrogen spectrum and the prediction of the antimatter [3] , and predated the experimental dis¬ 
covery of positron. In different parameter limits, the Dirac equation collapses to the Pauli equation, the 
Schrodinger equation, and the Weyl equation, respectively. Since the first experimental realization of 
graphene in 2003 ms\, much attention has been drawn to the study of the structures and/or dynamical 
properties of graphene and graphite as well as two dimensional (2D) materials [34], in which the Dirac 
equation plays an important role. This remarkable experimental advance renewed the interest on the 
theoretical analysis and numerical simulation of the Dirac equation and/or the (nonlinear) Schrodinger 
equation without/with external potentials, especially the honeycomb lattice potential [21 12011^ . 

After proper nondimensionlization and dimension reduction, the (linear) Dirac equation for the 
spin-1/2 particles with external electromagnetic potential in d (d = 3,2,1) dimensions reads |51I71ITB1 
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[niniiiniiMiisz]: 


(1.1) 






i=i 




1=1 


a 4 


«'(t,x), xG 


with the initial data 

( 1 . 2 ) 


= 0 ,x) = ^'o(x), 


X G 


where i = t is time, x = (a:i,... ,Xd)^ G is the spatial coordinate vector, dk = {k = 

l,...,d), 'll := 5'(t,x) = {tpi{t,x),tp 2 it,x),^ 3 {t,x),ip 4 {t,x))'^ G is the complex-valued vector wave 
function of the “spinorfield”, £ G (0,1] is a dimensionless parameter inversely proportional to the speed 
of light. In is the n x n identity matrix for n G N, Id := V{t,x.) is the real-valued electrical potential 
and A := A(t,x) = {Ai{t,x),... Ad{t,K))'^ is the real-valued magnetic potential vector. In addition, 
the 4x4 matrices ai, a 2 , and /3 are defined as 


(1.3) ai = 


0 (Tl 
(Tl 0 


a2 = 


0 (72 

<72 0 


as = 


0 (73 

(73 0 


/3 = 


0 —lo 


where ai, CT 2 , (73 are the Pauli matrices given by 
(1.4) 


(7l = 


0 1 
1 0 


(72 = 


0 -i 
i 0 


(73 = 


1 0 
0 -1 


The Dirac equation (HU) with (HU conserves the total mass 


(1.5) ||4((t,.)||2 := f |4r(t,x)pdx= / ^ |V' 2 (t,x)pdx= ||4'(0,-)f = ll^'o 

Jr'( Jr’>' 


t > 0 . 


In addition, if the external electromagnetic potentials are time independent, i.e. V (t, x) = D(x) and 
Aj{t,x) = Aj(x) {1 < j < d), the Dirac equation (II.II) with (11.21) conserves the energy 


(1.6) E{t) := 


/R'i 


-f D(x)|4'p - ^ Aj(x)4'*aj4' ] dx = ^ 1 ( 0 ), t > 0 , 

^1=1 ^ 1=1 


where denotes the conjugate transpose of In the nonrelativistic limit regime, i.e. 0 < e <C 1, 

as proven in the solution ih of the Dirac equation HU can be split into the 

electron part and the positron part, i.e.. 


(1.7) 



VC 




^2 

-G e**/"' 

0 


0 




Vo^ 


^4/ 


-f 0{e) = -f 0 (e), 


where both the ‘electron component’ and ‘positron component’ 4>p satisfy the (different) Schrodinger 
equation [ 121137 ) . In addition, a higher order 0 (£^) approximate model of the Pauli-type equation was 
provided in [ 3 T 1 | 32 ]. For details of the nonrelativistic limit of the Dirac equation HU, we refer the 
readers to [T 2 l[ 3 Tl[ 32 l[ 37 ] and references therein. 
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In practice, for lower dimensions d = 1,2, the Dirac equation (ED can be split into two equivalent 
sets of decoupled equations with two components each [7] and thus can be reduced to the following 
equation for $ := $(t,x) = x), x))"^ S as 


(1.8) i5t$(t,x) = 

with the initial data 


d ^ d 

+ ^^"3 + V (t, x )/2 - 
i=i 1=1 


$(t,x), xeR'^, 


(1.9) $(t = 0,x) = $o(x), xeR'^, d=l,2, 

where $ = (ipi, 1 ^ 4 )'^ (or $ = {ip 2 , in one dimension (ID) and under the transformation X 2 —t —X 2 
and A 2 —>■ —A 2 in 2D). As a result of its simplicity compared to (11.11) . the Dirac equation (11.81) has 
been widely used when considering the ID and 2D cases [71 15^1^ . 

Similarly, the Dirac equation (II.8p with (jl.9p conserves the total mass 

2 

(1.10) ||$(t,.)f := /" |$(i,x)pdx= / ^|(/)j(t,x)pdx= ||$(0,-)f = ||$of, i>0. 

Furthermore, if the external electromagnetic potentials are time-independent, i.e. D(<,x) = D(x) and 
Aj{t,x) = Aj(x) (1 < j < d), the Dirac equation (11.81) with (11.91) conserves the energy 


(1.11) E{t):= 




d ^ d 

---h l^(x)|$|^ - Y I dy: = E{ 0 ), t > 0. 

^ 1=1 ^ 3=1 


For the Dirac equation (11.81) . one can obtain the nonrelativistic limit which is similar to (ED and the 
detail is omitted here for brevity EI5TII571I57]. 

There have been extensive theoretical and numerical results for the Dirac equation (ED (or (ED) 
in the literatures. Along the analytical front, time independent states and dynamical properties have 
been thoroughly investigated, such as the bound states [18] , semi-classical limit [23] and nonrelativistic 
limit [1^331137] . etc. Along the numerical front, various finite difference time domain (FDTD) methods 
[UmillSIinnillS], time-spUtting Fourier pseudospectral (TSFP) method [TT1|^ and Gaussian beam 
method m have been proposed to solve the Dirac equation (ED (or ED)- However, most existing 
numerical methods are designed for the efficient and accurate simulation of the Dirac equation (ED 
(or ()1.8l) i in the parameter regime e = 0(1). In fact, for the Dirac equation in the nonrelativistic 
limit regime, i.e. 0 < £ <C 1, based on the theoretical analysis Eimi^isiHssiisziiin] , the solution 
exhibits propagating waves with wavelength 0(£^) and 0(1) in time and space, respectively. This rapid 
oscillation in time brings significant difficulties in designing and analyzing numerical methods for the 
Dirac equation (11.11) (or (11.81) 1 when 0 < £ <C 1. Recently, we have rigorously analyzed and compared 
the frequently used FDTD methods and TSFP method for the Dirac equation in the nonrelativistic 
limit regime [7] and shown that the meshing strategy (or £-resolution) for the FDTD methods and 
TSFP method should he h = 0(^/e), t = 0(£3) and h = 0(1), r = 0(£^), respectively, where h is the 
mesh size and r is the time step. Thus, the existing FDTD and TSFP methods are capable to solve 
the Dirac equation (ED (or (11.81) 1 efficiently and accurately in the regime £ = 0(1), and are much less 
efficient and time consuming in the nonrelativistic limit regime 0 < £ -C 1. 

The main aim of this paper is to design and analyze an efficient and accurate numerical method 
for the Dirac equation (ED (or (11.81) 1 which is uniformly accurate (UA) for e S (0,1]. The key 
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ingredients include adopting a multiscale decomposition of the solution of the Dirac equation m at 
each time interval with proper transmission conditions at different time intervals and applying the 
exponential wave integrator (EWI) with appropriate numerical quadratures which have been widely 
explored in solving highly oscillatory ordinary differential equations (ODEs) [JSIIMIEI] and dispersive 
partial differential equations (PDEs) [6l|9l[T0]. Then by a careful study of the error propagation and 
using the energy method, we establish two independent error estimates via two different mathematical 
approaches for the MTI-EP method as ^ and with mg depending on the regularity 

of the solution. Thus the MTI-EP method converges uniformly in space and time with respect to 
0 < e < 1. We remark here that a similar MTI-EP method has been recently designed and analyzed 
for the nonlinear Klein-Gordon equation in the nonrelativistic limit regime [S]. 

The rest of the paper is organized as follows. In section[2 we introduce a multiscale decomposition 
for the Dirac equation (11.811 and design the MTI-EP method. In section [H we establish rigorously 
error estimates for the MTI-EP method. Section |4] is devoted to the numerical results of the MTI- 
EP method and convergence rates of the solution of the Dirac equation to its limiting models in the 
nonrelativistic limit regime. Einally, some conclusions are drawn in section [SJ Throughout the paper, 
we adopt standard notations of Sobolev spaces and their norms, and use the notation p ^ qto represent 
that there exists a generic constant C > 0, which is independent of time step t, mesh size h and e, such 
that IpI < Cq. 

2. The MTI-FP method. For simplicity of notations, we shall only present our method and 
analysis for the Dirac equation dn in ID, while all the notations and results can be easily generalized 
to (11.81) in higher dimensions (2D) and (11.11) without any extra work. Denote 

(2.1) T =—ieaidx + cr^, W{t,x) = V{t,x)l 2 — Ai{t,x)ai, a; G M, 

where the domain of the operator T is and then the Dirac equation (11.81) in ID can be 

written as 

( 2 . 2 ) idt^{t,x) = 

We note that T is diagonalizable in the phase space (Fourier domain) and can be decomposed as 

(2.3) r=Vi- e^A n+ - - £2a n_. 


where A = d^x is the Laplace operator in ID and I is the identity operator, 11+ and n_ are projectors 
defined as 


(2.4) 


n+ = i 


2 1 


I2 + {I- e^A) 


- 1/2 


r 


n_ = - 

2 1 




Here \/I — e^A is understood in the Fourier space by the symbol \/T+~e^^ (^ G R) with domain (R). 
When A, I and I — e^A act on vector function $ = (^ 1 ,^ 2 )^, we mean that A, I and \/I — e^A 
act on two components (/)i and (/) 2 - It is easy to verify that n+ -|- n_ = I 2 and n+n_ = n_n+ = 0, 
n+ = n±. In addition, T and n± can be easily calculated in Fourier domain. 

As will be shown in the subsequent discussion, this decomposition of T is the key step for designing 
the uniformly accurate numerical scheme. The idea is to decouple the solutions into the eigenspaces of 
T in phase space by applying the projectors n+, and deal with the projections separately. It is worth 
noticing that the above formulation is in the whole space, while it is necessary to truncate the problem 
onto a bounded domain in computation. Thanks to the Fourier series, the operator T and its associated 
eigenvalue decomposition n+ can be explicitly computed in the phase space for the bounded domain 
case (see ()2.21L (12.221) '). which is crucial for the success of our method introduced below. 
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2.1. Multiscale decomposition. In order to design a uniformly accurate numerical method for 
the Dirac equation (11.81) (or (EB), from the experience in designing uniformly accurate methods for the 
nonlinear Klein-Gordon equation in the nonrelativistic limit regime [^ [T4lll9) . recalling that there exist 
propagating waves with O(e^) wavelength in time, a multiscale decomposition should possess O(e^) 
accuracy, so that the first order time derivative of the residue is bounded and a uniformly accurate 
scheme can be obtained. Thus, the first order Schrodinger decomposition dni) is inappropriate, and the 
second order Pauli-type decomposition (see [1^311133) 1 might work. However, due to the linearity of the 
Dirac equation (11.81) (or (HE)), we have a direct and better decomposition by applying the projectors 
n± to the Dirac equation HE HU- 

Choose the time step r := At > 0 and denote time steps as := nr for n > 0. From t = to 
tn+i, the solution $(t,x) = $(t„ -|- s,x) (denote $”(x) = $(t„,x)) to the Dirac equation (11.81) can be 
decomposed as 


(2.5) <i>(t„ -I- s,x) = e ^'I')j:"'(s,x) -I- 'I')_’"'(s,x)^ -I- ^4'^"(s,x) -I- 'I'i’"'(s,x)^ , 0 < s < t. 


where x), d(')j"(s, x)^ 


solves the coupled system for x G M and 0 < s < r as 


f x) = ^ (VZ-e^A - I) ^^'"(s, x) + n+ (lF^'i’^(s, x) + W^b"{s, x)) , 

(2.6) < x) = ^ (-VZ-e^A - I) ^')_’”(s, x) -k H. (lF«'Z"(s, x) -k lT«'l_’”(s, x)) , 

[«'+’"(0,x) = n+$"(x), «'l_’”(0,x) = 0, 


with W := W{tn -f s,x), and similarly ^d>^’"(s, x), d>^’"(s, x)^ satisfies 

f x) = ^ (VZ-e2A + I) ^^'"(s, x) + n+ (lF^'^”(s, x) + lF^'i’”(s, x)) , 

(2.7) < icls'I'^"(s, x) = ^ (—\// — e^A -I- /) ^/^’"(s, x) + n_ ^lT5'^"(s, x) -I- lF4'^’"(s, x)^ , 

[«'+’"(0,x) = 0, (0, x) =n_$’"(x). 


Following the analysis in [T^], it is easy to verify that 'k)|l"(s,x) = 0(1), 4'^’”(s, x) = 0(1),'I'l_’”(s, x) = 
0(£^), 4'^’"'(s, x) = 0(£^), and = 0(1) for A: = 1,2. Thus $(t„+i,x) can be evaluated numer¬ 

ically by solving the two coupled systems (EE) and (1^ via the exponential wave integrator Fourier 
pseudospectral (EWI-FP) method |6l[8] through the decomposition (12.51) . 


2.2. The MTI Fourier spectral discretization. As a common practice in the literatures HB 
I13II24II25II28II36II41II42) for practical computation, the Dirac equation (11.81) with d = 1 is usually truncated 
on a bounded interval H = (a, b) for $ := $(t, x) G C^, 


(2.8) idt^{t,x) = —T^{t,x) + W{t,x)^{t,x), x G H, t > 0, 

with periodic boundary conditions and initial condition as 


(2.9) $(t, x) is (6 — a) periodic in X, t > 0; $(0,x) = <i)o(x), x G H; 

where we use the same notations T and W{t, x) here as those in the whole space case (12.1|) by abuse of 
notation and we remark that the domain of There is (i7p(D))^ with Hp{Q) = {m G H^{n)\u{a) =u{b)}. 
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Then the systems (1^ and ism for the decomposition (12.51) with x G and 0 < s < r collapse to 


r x) = ^ {Vl-e^A - I) x) + n+ (w-^Yis, x) + lT^L’”(s, x)) , 

(2.10) i x) = ^ {-VI-e^ A - l) x) + n_ x) + x)) , 

is (6 — a) periodic in X, d>^’”(0,x) = n+$(t„,x), x) = 0, 

and 

f x) = ^ (Vi-e^A + I) x) + n+ x) + lT^^_’"(s, x)) , 

(2.11) i x) = j, (-VZ-e^A + I) «'!,’”(s, x) + n_ (lT«'^”(s, x) + lT^'!.’"(s, x)) , 

['I'^"(s,x) is (&—a) periodic in X, ^'^’"(0,x) = 0, ^'^’"(0,x) = n_$(t„,x). 

Choose the mesh size h := Ax = with M being a positive even integer and denote the grid 
points as Xj := a + jh for j = 0,1,..., M. Denote Xm = {U = (t/o, Di,..., Cm)^ | Uj G C^,j = 
0,1,..., M, Uq = Um} and the V norm in Xm is given by 

M-l 

( 2 . 12 ) \\U\\%=hJ2m^ UgXm. 

j=o 

Introduce 

Ym = Zm'xZm, with Zm = span|(;ii(x) = /r; = Z = - y,..., y - l|. 

Let [Cp{a,b)V be the function space consisting of all periodic vector function U{x) : [a, 6] —>■ C^. For 
any U{x) G [L^(a, &)]^ and U G Xm, define Pm ■ [L^(a, &)]^ —t Ym as the standard projection operator, 
Im ■ [Cp{a, 6)]^ —>■ Ym and Im '■ Xm —t Ym as the standard interpolation operator, i.e. 



M/2-1 

M/2-1 


(2.13) 

{PmU){x)= Y. 

{ImU){x)= Y 

, a < X < b. 


l=-M/2 

l=-M/2 


with 




(2.14) 

Ui = f U{x)e-^^‘^^-‘^Ux, 

0-a Ja 

^ M-l 

Ui= Y Uj 1 ^ 

j=0 

M M 

~~2' " ' ’ ”2 


where Uj = U{xj) when 17 is a function. The Parseval’s identity implies that 

(2.15) \\Im{U){-)\\l- = \\U\\i2, WgXm. 

The Fourier spectral discretization for (I2.10I) - (I2.11I) reads: 

Find := vl/^’;^^(s) = 'I/^’”^(s,x) G Ym (0 < s < r), i.e. 

M/2-1 _ 

(2.16) ^f^’”^(s,a;) = ^ a < x < b, s > 0, k = 1,2, 

l=-M/2 
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such that for a < X < 6 and 0 < s < r 


(2.17) 


and 


(2.18) 


^ (-VT^7^-I) d/i’"^(s) + n_ (vhvi/^’"^(s) + w^]:-^(s)), 

.«'+’"m(0) = Pm iTl+^tn,x)), 'Sh^MiO) = 0, 

^ {VT^^ + 1 ) vi/^-;^(s) + n+ + Vhvi/2_’"^(s)), 

^ (-vT^^ + /) 'I'!.’”m(s) + n_ + Vhvl/2_>"^(s)) , 


.«';’"m( 0 ) = 0, 'h!.’”M( 0 ) = Pm (n_<I>(t„,x)). 

We then obtain the equations for the Fonrier coefficients with 0 < s < r as 


(2.19) 


and 



zds (d/i”) ^ (s) = ^ 72 ('hi") /s) + n+ (iF^) i^) + n+ (iFvh^) ^ (S), 
*a.('hi’")^(s) =-^72('hi’”)^(s)+ nr(vFvhi3^)^(s)+ nr(iFvi/i"^)^(^ 

td.^){s) = -^7,(vi; 2_.")^(,) + nr(iFd/i3^)^(s) + nr(iFvhi;^^ 


(2.21) Si = + S+ = 5i + l, Sr=Si-l; 

and 11)*" and 11)" are the corresponding Fourier representations of the projectors n± as 

.V 


( 2 . 22 ) 


K = 


' l+Si 

2Si 

EfJ-l 


SAC 

25; 

2 ^2 
£ Mi 


/ _£i£i' 

nr = 2<5i(5i + l) 25, 

^ I _£i!± l+<?i 

25, 25, y 


\ 


2(5i 2(5i((5i + l)/ 

Using the variation-of-constant formula, we can write the solution as 


_ M M 


('hi")^(s) = e-'^r£/£^(vl;b")^(0) _i^*e-*^r(£-»)/£=n+ (^(iuvl/i;^)^(n;) + (lUvI/i^^)^(rc)^ dw, 


(^^•")^(s) = 


^7(5+ (s-u 


= g-i-S,"*" s/e 


— A / e 
Jo 


—iS"^ (s—w) /e^ 


)/^^nr (^(lFvlAi3^)^(u;) + (wYd%)^iw)j dw, 

t ((w^'hirM)/^) + (w^'hi-M)/^)) dw, 

(«'!.’") ^s) = e''^re/e"(^^")^(0) -ij^ g*5f ^ (^y) + (lF^'i”j^^) ^ (w)^ dw. 




’('hi")iO) 
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Using the initial conditions and choosing s = r, we can approximate the integrals via the Gautschi type 
quadrature rules [8l[22l|26l|23 or EWI [6HT0]. Using Taylor expansion and the equations (I2.17I) - (I2.18I) 
to determine the first order derivative, we can approximate the first integral in the above equation as 


(2.23) 


where 

(2.24) 


= Pi{T) n+(tuvb^-"^)^(0) + gr(r) n+ + , 


Pi (t) 


(&)• 


di {r) 




e 



sine 



and sinc(s) = with sinc(O) = 1. Note that (r) = 0(t) and (r) = 0{t^) for I = —. ■., ^ —1, 
and for the special case I = 0, Pq(t) = —it and q^ir) = —i\- Similarly, the other integrals can be 
approximated as 


(2.25) -l ^ ^ (zc)) dw 

^ nr(iu4^i"M)^(0) -^nr (^9,(tuvi/i"^)^(0) + a.(iuvi/^3^)^(0)) , 

(2.26) ^ (y;) + ^ (ic)) dw 

^p+(r) n+(fuvl/2j"^)^(0) + 9 +(r) n+ (9.(iU4/^-"^)^(0) + 5«(lUvI/!.’"^)^(0)) , 

(2.27) _ (-)/-= nr (^(fUvI/^"^)^(u;) + ^{w)^ dw 

^ nr /O) - ^ nr (a. ^ (o) + a« ^ (o)) , 


with c denoting the complex conjugate of c and 


(2.28) 


P'1 (r) = —ire 2 e^ sine 


^h. 

2^2 


dti^) = - 


re" 

'W 


I — e 2 e^ sine 


^h. 

2e2 


For simplicity of notations, by omitting the spatial x variable and denoting 

(2.29) /l’"(s) = fU(t„ + s)vl/^’"^(s), /l’"(s) = 5l’”(s) = d.Wit^ + s)'I/^’;)^,(s), 

in order to design a uniform convergent method with respect to e G (0,1], we find the solutions should 














be updated in the order from small component to large component as 


I ^ -P^Wnr(/i’")/0) - 9 ,+ (r)nr(gi”)^(0) - g+(T)nr ((/i’”)/0) + (/!’”)/O)) , 

|(vI/^+”)^(r) « p+(r)n+(/!■") ^(0)+g+(r)n+(5!.’") ^(0) + g+(^^ ((/^’")^(0) + (/!-'^)^(0)) , 

with initial values and derivatives determined from (I2.17p - (l2.18p as 


('tt”)/0) = n+( 5 (tj)„ ('i'i-”)^(0) = 0, (i'i")^(0) = 0, ('I't")^(0) = nr( 5 (tb)„ 

(a®i’”„)/o) = -i n,-(w((„)»i'"„(0))_, (a,»i"„)^(0) = -i n+(w'(t„)<ti'j,(0))^, 

^ -. nr(tr«„l^„(0))^, 1 = -f.f -1. 


where the derivatives i9s'I'+’m( 0) and are approximated using filters 2 sin(/r^T/2)/T instead 

oi fif {I = ^ — 1) to avoid loss of accuracy, then followed by 


(vE-^-")^(r) « +p-(r) n+(4-")^(0) + n+(gi’")^(0) 

+qr{r) n+(^(/i’-)^(o)+ (/!’-)^(0)), 

(vE-2j")^(r) _^nr(/!’”)^(0) -^nr(5-’”) 

nr ((/+”)/O) + (/-”)/0)) ’ 


with (0) and (0) approximated in another way as 




2.3. The MTI-FP method in ID. In practice, the integrals for computing the Fourier transform 
coefficients in (12.141) . (I2.23|) - (|2.27p are usually approximated by the numerical quadratures. Let be 
the numerical approximation of the exact solution to the Dirac equation (12.81) for n > 0 and 

j = 0,1,..., M, and denote $” £ Xm as the numerical solution vector at time t = tn, in addition, let 
^fc,ra+i 2 ) i^e ^;]^e numerical approximation of 4'^’"(r, Xj) tor j = 0,1,..., M and n > 0, and 

denote f/” = V{tn,Xj), A^j = Ai(tn,Xj) and W" = W(tn,Xj) = — ^ij^i for j = 0 , 1 ,..., M 

and n > 0. Choosing <I>° = <I>o(a;j) for j = 0,1,... ,M, then the multiscale time integrator Fourier 
pseudospectral (MTI-FP) method for discretizing the Dirac equation (11.81) in ID reads for n > 0 and 
j = 0,1 ,..., M as: 


(2.30) $ 


n+1 


,-iT/e^ 


4- 


l.n+l 
+ J 


4- 


L,n+l\ 


r/e^ 


^4/^"+l q- 


M/2-1 

^ ($"+i),e 

i=-M/2 
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where 


(2.31) 




fc,n+l 


MI2-1 

E 

l=-M/2 






k = 1,2, 


with 

(2.32) 

= -EEnrCm, -EEnrSo, -EEnr , 

+ 5 +(^)n+(^^ + Q+(r)n+ , 

(virh-+i)^ =e-^'^^+p-(^)n+(7r5^ + g-(^)n+(^^ + g- (^[fl)^ +[fb*)^ , 

^ g*-^ _p- (^)n- (/2)^ - q- (r)n-(^2 )^ _ q- (r)n- ^(/+*)^ + (/-)^^ , 


and 

(2.33) 

M/2-1 

/±., = E 

1^-MI2 
M/2-1 


M/2-1 


il,* 


/^— iyi/^ — ± - --—-- 

f±,j= E /n 5l, = E 

i=-M/2 ^ 2 i i=-M/2 

M/2-1 


! •^ — -*-/ \ ±v±! ^ — j-/ 

/-.; = E (E’*) //*. = ^ (pA^^m{x,-a)^ j = 0,l,...,M, A = 1,2 

' _A/f /9 1 = — Mj 2 ^ 


Z=-M/2 


with 

(2.34) 


^, = n+(i^„ (4'L), = 0, (n)i^ ('f'l), = nr(i^i, 

(n)^n+(4),, (^l)^ = -*nr(E),, / = - 

('i'lX = - *nr(7/j,, 

fl^ = 1174/1,^-, fl^ = ffl, = ftll^(t„, x,) 4 /|_^-. 


M 
2 ’ 


/li= 


7 M'" (*-7‘) 


n2,* _ 

J+,3 ~ 




M 


- 1 , 


j = 0,l,...,M, fc = l,2. 


Remark 2.1. When the electro-magnetic potentials in the Dirac equation (12.Sp are time indepen¬ 
dent, i.e. V{t,x) = V{x) and Ai{t,x) = 74i(a:), the above scheme (12.30p - (12.341) will be simplified with 
g± j = 0, k = 1,2, j = 0,1,..., M, i.e., the {g±)i term vanishes. 

Note that in the above MTI-FP method, we first evaluate the small components and 4'/"”'’^, 

then approximate their time derivatives via finite difference approximations, and finally use them in 
the evaluation of the large components and 

This MTI-FP method for the Dirac equation (II.Sp (or (11.11) 1 is explicit, accurate, easy to implement 
and very efficient due to the discrete fast Fourier transform. The memory cost is 0{M) and the 
computational cost per time step is 0{M log M). As will be shown in the next section, it is uniformly 
convergent in space and time with respect to e G (0,1]. 
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3. A uniform error bound. In this section, we establish two independent error estimates for 
the MTI-FP method (12.301) via two different mathematical approaches. Let 0 < T < oo be any fixed 
time. Motivated by the results for the Dirac equation (jl.8l) (or (HU) in usm, we make the following 
assumptions on the electromagnetic potentials 


(^) II^IIiV2.°°([0,T];(1Vp'"“'“)2) + II^i|Im/2,cx=([o_T];(w;"0'°°)2) ^1) nio > 4, 

and the exact solution $ := of the Dirac equation (12.81) with e G (0,1] 




5, L 




< — 

^ £ 4 ’ 


where iL™(D) = {u \ u G dl.u{a) = dh^uib), I = 0,... ,m — 1} and Wp’°°{il) = {u \ u G 

dlu{a) = dlu{b), I = 0,..., m — 1} for m G N. We remark here that the assumption (B) is 
equivalent to the initial value $o(x) G [T^[55] under the assumption (A). 

Theorem 3.1. Let $" G Xm be the numerieal approximation obtained from the MTI-FP method 
p.30[l and denote ^f{x) = Im(^"'){x) G Ym- Under the assumptions (A) and (B), there exist constants 
0 < To, ho < 1 sufficiently small and independent of e, such that for any 0 < e < 1, when 0 < r < tq 
and 0 < h < ho, we have 


(3.1) ||$(t„,.)-$ 7 (.)||^. + 0<n<-, 

T 

which yields the uniform error bound by taking minimum among the two error bounds for s G (0,1] 

(3.2) \mtr,,-)-<^fi-)\\L 2 <h'^«+ min \^,T^+eA<h^°+T, 0 < n <-. 

o<E<i j T 


Remark 3.1. From the analysis point of view, we remark that the assumption in (A) is 

necessary such that the exaet solution ^{t,x) of the Dirac equation remains in which would 

give the spectral accuracy in space. In practice, as long as the solution of the Dirac equation (or 

m is well localized such that the error from the periodic truncation of potential term W{t, x)^{t, x) 
is negligible, the error estimates in the above theorem still hold. 

In the rest of this paper, we will write the exact solution $(<, x) as $(t) for simplicity of notations. 

M/2-1 __ 

Define the error function e"(x) = ^ g Ym for n > 0 as 

i=-M/2 

(3.3) e’^(x)=PM($(in))(cc)-$?(x)=PM($(t„))(a:)-/M($”)(x), x G D. 

Using the assumption (B), triangle inequality and standard Fourier projection properties, we find 

(3.4) ||<i>(t„, •) - $7(.)|U2 < ||$(t„,.) - PM(<i>(i„))(-)llL2 + ||PM($(i„))(-) - $?(-)IIl2 

</i™“ + ||e"(-)||L 2 , 0<n<-. 

T 

Hence, we only need estimate ||e"'(-)||i 2 . To this purpose, local truncation error will be studied as the 
first step. Since the MTI-FP method (12.301) is designed by the multiscale decomposition, the following 
properties of the decomposition (I2.10I) - (I2.11I) are essential for the error analysis. 
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From t = tn to tn+i, let x) (0 < s < t, k = 1,2) he the solutions of the systems (12.101) and 

(12.111) . and the decomposition (12. 5|) holds as 

(3.5) + s,:c) = x) + x)'^ + e*"/"' x) + cr)) , x G ft. 

Then the error e"+^(x) (n > 0) (13.31) can be decomposed as 

(3.6) e"'''^(x) = ^z)^"''’^(x) + zL’"''’^(x)^ + ^z^’”''’^(x) + z^’"'''^(x)^ , x G 11, 


with 


m/2-1 

(3.7) z^’"+i(x)= ^ (z^’"+i);e*^*("-“)=PM(^±’"(T))(x)-JM(^±’”+')(:r), fc = l,2. 

J=-M/2 


By the same arguments in [T^], we can establish the regularity results. 

Lemma 3.2. Under the assumptions (A) and (B), the exact solutions dt^’”(s,x) (0 < s < t, 
0 < n < ^ — 1) with k = 1,2 of the systems (12.101) and (12.111) satisfy 


(3.8) 

(3.9) 

(3.10) 


ll^±”llL~([0.r];(ff™O)2) 

ll^+”llL“([0,r];(i/;*i)2) 


< 1, ||3ss'1' + ”||loo([0_^];(//™4)2) + |||9ss'I'^’"||loo([0_t-];(//™4)2) < 1, 
)2) ^ ll^s^i^+"llL“([0,r];(ff“'‘)2) + II^ss^'1i"||loo([o,t-];(//^4)2) < 

+ ||'I'i’”|lLoo([o_T-];(//^i) 2 ) < mi = mo-I, / = 1,2,4. 


Proof Noticing the properties of the projectors n± and the assumption (B), we get that the initial 
data 'P!(.(0,x) G 9s'P!(.(0,x) G with uniform bounds. The estimates for and 

9shave been derived in m, where one only needs to replace the whole space Fourier transform 
with the Fourier series on torus. Thus, the proof is omitted here for brevity. It remains to estimate 
9ss4)'^’"'. Here, we show the case fc = 1, while the case fc = 2 is quite similar and the detail is omitted 
here for brevity. Differentiating (12.101) with respect to s, we obtain for 9ss'I'£"'(s), 

*9ss^'i”(s) =- (/-e"A)-^/' Aas^i”(s)+n+as (vF^i”(s,x) + VF^^”(s,x)) , 

i9ss^'i’"'(s) =- ^ 2 ^ —^s^-"(s) + n+9s (^VF4'+"(s,x) + VF5'l_’"(s,x)^ . 

Since for any $ G with 2 < m G N, we have 


(3.11) 


(l-e^A) ^^^A$ 


HI 


,-2 < 


Hff! 


y/I - e2A + / 


$ 




H” 


which implies the bounds (13.91) for 9ss'l'^’”(s), in view of the estimates for and assumption (A). □ 
Having Lemma 13.21 and the decomposition (13.6|) . we can define the local truncation error ^^’”(x) = 
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M/2-1 

E 

l=-M/2 

(3.12) 


(x G 11, A; = 1,2, n > 0) for the MTI-FP method (|2.30l) - l|2.34l) as 


+ g+(r)nr(^), + g+(r)nr M/KO))^ + (/1(0))^ 

= (vl/f (r))^ -p+(r)n+(H(5), - g,+ (r)n+(^), - g+(r)n+ N/^(O))^ + (/EO))^ 

(C+")j = («'+"(r))^ - -p-(r)n+ (/i(0)); - {t)U+ (gi(0))/ 

/ — _M M _ 1 

L 2 1 1 9 -*-5 



/ii 


-9/ ('r)n/' 


(/-’*( 




+ Pz (EH; (/^(0)); + g; (EH; (ff2(0)); 


+9/ (^)n. 


/+’*(-) +(/^ 



where 


(vi/),--)^(o) = n+($^),, (4 -E)^(o) = o, (vi/2^”)^(o) = o, (vi/E)^(o) = n-(¥^)„ 
(3-13) <; /|(s) = lF(t„ + s)d/|’"(s), /|(s) = lP(t„)vi/|’”(s), g|(s) = 

/l'*(s) = iTP(t„) (vI/E+'(s)) , /+’*(s) = iVF(t„)(vI/^’"+'(s)), 0<s<T, A = 1,2, 


and 


M/2 __ 

(3.14) ^|’"(s,x)= ^ 

l=-M/2 


pi/ii(x-a) 


, X G n, 0 < s < T, A = 1,2, 


with 


(3.15) 


(ipW), 


_ ,-2sin(/ifr/2) 






(vi/i:"(s))^ - zn+(;^),, (4 /E(s))^ = -^^-(f^))^ 


(s))^ = -*n+(/2(s))^, (vi/!.’"(s)) 


„-2sin(;i,^T/2) 


I * ^ 


(vi/2j"(s))^-*nr(H5),. 


We have the following estimates for the above local truncation error. 

Lemma 3.3. Under the assumptions (A) and (B), there exist constants 0 < ro,Ao < 1 sujficiently 
small and independent of e, such that for any 0 < e < 1, when 0 < r < rg and 0 < h < ho, we have the 
error estimates for the local truncation error G Ym in (13.121) 

(3.16) ||^^■"(.)||^. + ||^fc■"(.)|l^. <r(A™»+r2+£2), 0 < n < ^, A = 1,2. 

Proof. We will only prove the estimates p.l6l) for A = 1, as the proof for A = 2 is the same. Using 
the fact > 1 and the definitions of p^{t) and qfir) {I = —M/2,... , M/2 — 1), we have 


(3.17) 


\pH'^)\ - U ^ \pti'r)\ < e^, \q+{T)\ < TE^. 
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Multiplying both sides of the equations in the system (12.101) by e and integrating over fl, we 

easily recover the equations for which are exactly the same as (I2.19I) - (I2.20I) with being 

replaced by 

Following the derivation of the MTI-FP method, it is easy to find that the local truncation error 
comes from the approximations in the integrals (12.231) . (j2.25l) . (I2.26P and (I2.27p . In particular, for 
I = —M/2,..., M/2 — 1, we have 

(3.18) + ds+^)ur{f^)), 

+ dt ('^) nr {9lio))i + qt (t) nr 

(3.19) + ds-p-ir)U/-{f^))^ 

- qri'^) n,+ (5+(o)); - gr(^) n+ ^(/i(o))^ + (/l’*(r))^^ . 



Type I estimate. Here we prove the first kind estimate in (I3.16F Using Taylor expansion, we have 

f>T pS pSi 


(3.20) (e"),=-z 

(3.21) 




/O ^0 ^0 

nT pS pSi 


0 ^0 ^0 


^ ('9s2S2 f+{s2))i + {ds2S2f-{s2))i'j ds^dsids + (r?!);, 
(T-s)/e"jj+ {^[ds.2S2fl{s2)) 1 + {ds2S2f-{s2))^ ds2dsids + 


Mj2-l _ 

where q±{x) = ^ with 

l=-M/2 


=p+(r) nr (^-(/i(O)), + (4(0)),j + 4(r) nr (^-(4(0)); + (4(0)), 

-U 


+ 4(r) nr - /4 


(4), =Pi (r) n+ u/i(o)),-(4(0)),] + Qi (r) n+ ((4(o)),-(4(o)), j 



+ di (^) n( 


+ I (jhn 


,-(/i(o)),+(rr^)),-(/fc 


Here 4’"(s) is given in (12.2911 with ^'^’4 being replaced by Since ||n ,^||;2 < 1 (/ = —■ ■ ■, ^ —1) 
with IIQIIp being the standard P norm of the matrix Q, using (13.1711 and triangle inequality, we obtain 


l(4),l 


+ T^ 


+ T 


(4(0)),-(4(0)), 
(440))^-(4(0)) 


(4(0)),- (4(0)), 


+ T^ 


(/rr)),-(/uo)), 


+ r^ 


(440)),-(440))^ +4 (440))^-(4(0))^ 
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From the Parseval’s theorem, we get 


\W-{-)\\h - iMifimwh + A\PM{gXm - iMiaXmwh 

+ A\PM{fX^m - + A\PM{jX^m - 

( 3 . 22 ) +rX\iMux^m-iMUX{mA+A\iM{fX^m-iMUXm\^^^^ 

Recalling assumptions (A) and (B) and noticing Lemma 13.21 we have 

/i(0) = w{tn)A"io) e 9±{0) = dtWitr.)AA{0) e 

/i’"(0) = G 

Employing (j2.15l) and Cauchy inequality, for toq > 4, we can bound ||ry)_(.)||i ,2 from p.22l) as 




M-l 


.h-YX |hF(t„,a;j)(9s«'+"(0,a;j) - 
\ j=o 


M-l 


+ r" 


\f^Y \w{tn,Xj){dsA''iO,Xj) - A''{0,Xj)) 

\ i=0 

<r(/i™« + A + (|1/m(9.vI/^"( 0)) - 4'i’"(0)|U2 + WlMidsA^A) - 'hi’"(0)|U2) 

(3.23) <r(h™« + A + (||Pm( 5«4'^'"(0)) - ^+-"(0)|U2 + ||PM(9«^i’"(0)) - ^L’"(0)|U2) . 


Using the equation (12.171) . we get 


2sin(/ifT/2) 


{dsAlA) - AXnxo) = - *—(^('i^r(o)), - in^A), 


-in? 


-ifX 


. I 2sin(^fT/2) 


(5^"),(0) - (i?),(0) = -iA [CfXA)i - ifXAX) > 

and 

(3.24) ||PM(a.4/i’"(0)) - /M('hL’"(0))|U= < WPMifXA) - /m(4(0))|U2 < 

Noticing that | sin(s) — s| < ^ (s G K), we have 


4-. 


2sin(^fT/2) 




A 


1 2 sm{AT/2) 

2^' ; 


< 4 , M M 

/ = -Y,..., y-l, 


which leads to 


(a.vi/^’^),(o)-(vi/^'‘)(o) 


<- 


(vhi”(o)), - aXAo))i 


(/i(0)),-(/i(0)), 


+ Tai 


AXAo))i 
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and for toq > 4, 

||PM(a«vi/^-"(o))-/M(4^i"(o))|U^ 

< i||PMK’”(0)) - /M(4'i”(0))|U. + r||PMK’"(0))||ff4 + WPMifim - /m(4(0))|U2 

+T + h^° It. 

Combining the above estimates with (I3.23L we derive 

(3.25) hU-)\\L^ ^ + T^) + r2(h’"'> + /r + r) < + r^). 

By the same procedure, ||??+(-)IU 2 can be bounded as 

\\vU-)\\l^ <r{r^ + + T^PM{ds^]l^m - PM('hi’”(r))/r|U. 

+ r||PM('hi’”(T))-/M('h^”(r))|U2, 


where Taylor expansion gives 

a^vl/(j"(0) - vl/(j"(r)/r =-T / / dss'^-'^isiT)dsids. 

Jo Jo 

Thus, recalling Lemma l3.2[ we estimate 


(3.26) hl{-)\\L^ < r(r^ + h™«) + ^ (^h^o + . 

Now, Lemma EH] together with (13.201) . p.21|) . (I3.22p . (13.231) and (I3.26|) implies 

II^±"'(-)IIl2 <T^ (^\\dssiW{tn + s)4'+’”(s))||i=o([o_^].(i2)2) + \\dss{W{tn + s)4'(_’”(s)) ||ioo([o,rl;(L 2 ) 2 ) j 

(3.27) +||^i(.)||^, + 


Type II estimate. Now we prove the second estimate for ^^’"(a;) in (13.161) . Starting from ()3.18|1 and 
(13.191) . we treat the terms involving /:J_(s), /+(s) and g+(s) in the same way as in proving (13.271) . and 
leave the rest terms as 


SB, =(ci), - * r 

Jo 

SB, =(B),-* r 

Jo 


^utir-syp nr(/IB),rfs + g+(r) nr(/l(0))^, 
n+ (Hw) ^ ds - qr (r) n+ (/!’*(r)) ^, 


with C±(2;) 


M/2-1 

E 


,=-M/2 


(C±);e*w(^-“) 


satisfying 


llCi(-)llL^<r(h-'>+r2). 


16 























The proof of the above decomposition and the corresponding error bounds for (x) is identical to the 
proof of (13.271) , and thus is omitted here for brevity. Applying triangle inequality and (I3.17|) , we have 



< 

(C), 

^0 


ds + re^ 

(/KO)), 


< 

(O'), 


(m, 

ds + 

(/^)), 


Recalling Lemma [3^ which implies ^ we know ||/1(s)||l 2 < e^, ||/ 

It and 

ll/-(0)llff"‘O-l ^ ll^d (O)lli/'"0-l ^ l|-fM(/+(0))||//m(,-l < II/+(0)11 //mo < 1. 
Hence, using the Parseval’s theorem, we get 






<llri 


L2 


< 


L2 


L2 




+ T 


+ T fL 


L“([0.t];(L2)2) 


L“([0,r];(L2)2) 


+ Te^ 


+ T^ 


/m(/1(0)) <r{h"^o+T^+e^), 


L2 


/M(/i’*(r)) <r(h-'>+r2+£2), 


L2 


which, together with p.27p . completes the proof for (13.161) . □ 

Subtracting (12.321) from (13.121) . noticing (12.331) and (I3.13p . we get error equations for z^’”^^(a;) 
(fc = 1,2) in dSIl) as 


(3.28) 


(3.29) 


J (zi-"^), 

= + (d’” 


1(Z^'"+^), 

= {Tir\F\tf 


M/2-1 



'x)= E 

(/c = 1,2) is g 

i=-M/2 


(^-“), = 

-p+(r)nr(d’" 

-9/('r)nr 

II 

}E 

p+(^)n^(d’")^ 

+ g+(dn^(' 

(^), = 

pr(r)n+(A/”)^ 

+ dr ('^)ni+ (' 

(^-“), = 

-pr(r)nr(A!’" 

)j-dr(^)nr 


+ (Ci”) 

-2,n\ I //-2,n\ 


M/2-1 M/2-1 

with ^ gYm (k = 1,2), G±±’”(a;) = ^ ^ 

l=-M/2 l=-M/2 

F±±’”(x) = f; (yfe±.")^gzw(x-a) g Ym, and = f) (F^=F,*)^g*wG-a) g Ym {k+ = 

l=-M/2 l=-M/2 
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1, k- = 2) defined as 


(3 30 ) = {F^ni = {fm)i-if±)v {Gib'll = {9i^ mi 


V-+ )l 
~~k,n / 


For the electromagnetic error part (fc = l,2,0<n< ^ — l),we have the lemma below. 

Lemma 3.4. Under the assumptions (A) and (B), the electromagnetic error part F^"'{x) G Ym 
( k = 1,2, 0<n< — — 1) defined in (13.291) with p.30l) satisfies the bounds as 


,n+l / 


|F^±’”(-)||l= + + ||e"(-)|U 2 , k+ = 1, k. = 2, 

which implies that 

(3.31) II J-^’”(-)I|l^ < r(h™° + ||z^±’”+^(.)|U. + ||e"(-)|U^), || J-^’”(-)IIl 2 < r(h-“ + ||e"(.)|U.). 


Proof. Recalling the assumptions (A) and (B), Lemma [3^ (I3.30L (13.131) . (13.151) . (12.331) and (12.341) . 
applying the Parseval’s theorem, we have 


M-l 


\\Fl’^i-)\\h <\\lMiflm - iMifDWh = h^ \w{t„,x,)mf^i0,x,) - 4/^,,.) 

1=0 

M-l 2 

<h + II^M($(in)) - /M($")|li. < + ||e'^(-)|li., 

1=0 

and similarly we have ||A^’”(-)||l 2 < h'^° + ||e"(-)||i 2 . Using the same idea, we can obtain 

||Gi”(-)|U= + ||G^J”(-)IIl= <+ l|e”(-)llL^ 

and 

l|i^^*(-)llL= <l|/M(/l’*(r)) - /m(/-’*)||l 2 < -||/M('hL’"(r)) - 

T 

T 

ii4’*(-)iil= <ii/M(/i’*(r)) - + iiz^"+^(-)iil0- 

T 

It remains to estimate Again, in the same spirit of the above arguments, we arrive at 

(3.32) l|i^^’”(-)IU= < l|/M('h|’”(0)) - /M('hi)IU- 

Comparing (13.151) with (12.341) . noticing the properties of 6^ and the arguments in the above proof, v 
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find 


I1/m('&L’"(0)) - <I|/m(4(0)) - iMifDh^ < + ||e"(.)|U=, 

||/M('i'^”(0)) - <I|/m(/!(0)) - /m(/!)I1l^ < + ||e'^(-)IU^, 

||/M('i'i"(0)) - <^l|/M('I'i"(0)) - /M(«'i)||L= + I|/m(4(0)) - iMifDh^ 

<i(/i-o + ||e-(.)||^,) + /,-o + ||e"(.)||^,, 
r 

||/m(«''_’"( 0)) - /M('i'"_)||L^ <-I1 /m('p'j"( 0)) - /M('Pi)||L2 + \\lMif-iO)) - /m(/!)||l2 

r 

<-(/i™" + l|e"(-)l|LO + /i™" + l|e"(-)llL^, 

r 

which implies the bounds for ||F£’”(-)||i 2 in view of (|3.32p . Combining all the above results, recalling 
(13.301) and properties of the coefficients pfir) and q^ir) in p.l7p . we conclude that (13.311) holds. 0 
Now, we are ready to prove the main theorem. 

Proof of Theorem 1,9.11 Recalling the decomposition (13.61) and the error equation (I3.28|l , we get 




(3.33) 





Mjl-X __ 

with x” 

■(x) = 



i 

'.=-M/2 

(3.34) 

(^), 

^ g*T(2fe-3)/E^ 


fc=l,2 


M 

T’ 



In particular ||e°(-)||i 2 = ||Pm($o) - Im{^o)\\l^ ^ ■ 

Taking the P norm of the vectors in the error equation (13.281) . then summing together for I = 
—M/2,..., M/2 — 1, utilizing Lemma [3.41 and Parserval’s theorem, there holds for 0 < r < 1, 

<ll-^i’"(-)llL^ + lld’”(-)llL^ < + l|e”(-)llL0 + lld’"(-)llL^, 

Il4’”^'(-)IU= <II-^+’”(-)IIl^ + lie+’”(-)llL= < rih'-- + l|e”(-)llL0 + lie+’”(-)llL^, 

and so 

IIX”(-)II + l|e”(-)llL^ + + ||z^”+^(-)||l=) + E + lld’”(-)llL^) 

/c^l,2 


From Lemma 13.31 on the local truncation error we get 

(3.35) ||^-(.)||^, <^||e-(.)||^.+r(/r™«+TV£ 2 )^ 0 < n < - - 1, 

r 

(3.36) ||^»(.)||^, <^||e"(.)||^.+r(/r'"»+r2+£2), q < n < - - 1. 

T 
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Now, taking the P norm of the vectors on both sides of (13.331) . making use of the orthogonal properties 
of 11;^ where + |nj”vp = |vp for all v G 0 ^, 01,02 € M, we can have 


(e"+i), = n+(e-), + nr(e"), + |(x")/ 

+ 2Re ((e-*^‘/^^n+(^), + 


where Re{c) denotes the real part of the complex number c. Applying Cauchy inequality, we find 


(3.37) 


VA ^ ^ ^ /—\ 1 9 1 1 /—\ 9 Af M 

(e”+i)i - (e’^), <r|(e-);|2 + _|(^n)j2^ ^ i. 


Summing (13.371) together for / = — ..., ^ — 1 and using Parseval’s theorem, we obtain 


2 -+ 1 (.)|| 2 , _ ||e-(.)|||, < r||e"(.)||i. + -||x”(-)lli=, Q < n < - - 1. 

T T 


(3.38) 

I 7 

Summing (13.381) for indices 0,1,..., n and using (I3.35I) - (I3.36I) . we derive that for 0 < n < — 1, 

(3.39) \\e-+\-)\\h - l|e°(-)lli^ E N™(-)lli^ + nr 

n 

(3.40) \\e-+\-)\\h - I|e°(-)lli2 <r E \\e^{-)\\h + nr{h^'> + . 

m—1 

Since l|e°(-)llL^ < , Gronwall’s inequality will lead to the conclusion when 0 < r < tq < 1 and 

0 < h < ho < 1 sufficiently small 


(3.41) ||e-+i(.)||2, < ^ ||e"+i(.)||2,^ < (/i™o+r2+e2)2, q < n < ^ - 1. 

In view of (13.4L we conclude that (13.11) holds. □ 

4. Numerical results. In this section, we present numerical tests on our MTI-FP method (12.301) 
and apply it to study numerically the convergence of the Dirac equation (11.81) to its limiting Schrodinger 
model 113 and the second order limiting Pauli-type equation model. To this purpose, we choose the 
electromagnetic potentials in the Dirac equation (11.81) with d = 1 as 

(4.1) Ai{t,x) = ^2 ’ ^ft. = 1 ^, ^2 ’ a;GK, t > 0, 

1 -I- 1 -I- 

and the initial data $o(x) = ((/)i(x), ^ 2 ( 2 ;))^ as 

(()i(x) = e~^ P, 4'2 {x) = P, X G K. 


(4.2) 
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Table 4.1 

Spatial error analysis of the MTI-FP method for the Dirac equation in ID. 


eh.r(2.0) 

ho — 2 

hol2 

hol2^ 

hol2^ 

hol2^ 

£0 = 1 

1.65 

5.74E-1 

7.08E-2 

7.00E-5 

8.53E-10 

£o/2 

1.39 

3.45E-1 

7.06E-3 

6.67E-6 

9.71E-10 

eo/22 

1.18 

1.67E-1 

1.71E-3 

1.43E-6 

l.lOE-9 

eo/2^ 

1.13 

1.46E-1 

1.03E-3 

6.77E-7 

9.16E-10 

eo/2^ 

1.15 

1.45E-1 

8.52E-4 

4.86E-7 

1.33E-9 


4.1. Accuracy test. The Dirac equation (11.81) with d = 1, 6ID and (14.21) is solved numerically on 
an interval n = (—16,16), i.e. a = —16 and b = 16, with periodic boundary conditions. The ‘reference’ 
solution ^{t,x) = {(j)i{t,x),(j) 2 {t,x))'^ is obtained numerically by using the TSFP method [7] with a 
very fine mesh size and a small time step, e.g. he = 1/32 and Tg = 10“^. Denote 4))) ^ as the numerical 
solution obtained by the MTI-FP method with mesh size h and time step r. In order to quantify the 
convergence, we introduce 


ehAtn) = ||$"-$(t„,-)l|p 


\ 


M-1 

'■El*. 

3=0 


n 




Tab. 14.II displays the spatial errors Ch.rit = 2.0) with r = 10“"^ for different e and h; and Tab. 14.21 
lists the temporal errors eh.rit = 2.0) with h = 1/32 for different e and r. From Tabs. 14.1114.21 and 
additional numerical results not shown here for brevity, we can draw the following conclusions for the 
MTI-FP method: 

(i) For the spatial discretization error, the MTI-FP method is uniformly spectral accurate for all 
e G (0,1] (cf. Tab. 14.11) . 

(ii) For the temporal discretization error, the MTI-FP method is uniformly convergent with linear 

rate at 0{t) for e G (0,1]. For any fixed 0 < e < 1, when time step r is small, i.e. t < (upper 
triangle part of Tab. Ol) . second order convergence at is achieved; when e is small, i.e. e < t 

(lower triangle part of Tab. 14.2p . again second order convergence at O(t^) is achieved. However, 
near the diagonal part where r ~ (cf. the underlined diagonal part of Tab. 14.21) . degeneracy of 
the convergence rate and the uniform linear convergence rate for the temporal error are observed. In 
particular, the underlined errors in Tab. 14.21 degenerate in the parameter regime r ^ e^, which has 
been predicted by our error estimates (EH). 

4.2. Convergence of the Dirac equation (11.81) to its limiting models. Similarly to (HH, 
when e —>■ O’*', the solution of the Dirac equation (11.81) satisfies, i.e. first order limiting model 

(4.3) $(t,a:) = -|-e*‘/^Vpe 2 + 0(e), ei = (l,0)'^, e 2 = (0,l)'^, 

where c/g := </e(t,a') G C and (j)p := G C satisfy the Schrodinger equations [T^IMUTTIHU] . 


(4.4) 


idt4>e = 




idt4>p = 


-A + V{t,x) 




X € 


t > 0, 


and the initial data is determined through (EH. 
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Table 4.2 

Temporal error analysis of the MTI-FP method for the Dirac equation in ID. The convergence order is calculated 
as log 2 (e/i, 2 T/eh,r)- 



eh,.(2.0) 

To = 0.1 

to/2 

To/r 

to/2" 

to/2" 

to/2" 

to/2" 

to/2' 

to/2« 

£0 = 1 

3.69E-2 

9.18E-3 

2.29E-3 

5.73E-4 

1.43E-4 

3.58E-5 

8.94E-6 

2.24E-6 

5.59E-7 

order 

- 

2.01 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

eo/2 

5.98E-2 

1.51E-2 

3.77E-3 

9.45E-4 

2.36E-4 

5.90E-5 

1.48E-5 

3.69E-6 

9.23E-7 

order 

- 

1.99 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

2.00 

eo/2^ 

1.91E-1 

5.67E-2 

1.47E-2 

3.74E-3 

9.39E-4 

2.35E-4 

5.87E-5 

1.47E-5 

3.67E-6 

order 

- 

1.75 

1.95 

1.97 

1.99 

2.00 

2.00 

2.00 

2.00 

eo/2^ 

7.12E-2 

7.17E-2 

4.90E-2 

1.48E-2 

3.89E-3 

9.84E-4 

2.47E-4 

6.17E-5 

1.54E-5 

order 

- 

-0.01 

0.55 

1.73 

1.93 

1.98 

1.99 

2.00 

2.00 

eo/2"‘ 

1.78E-2 

1.76E-2 

1.80E-2 

1.82E-2 

1.22E-2 

3.73E-3 

9.79E-4 

2.48E-4 

6.21E-5 

order 

- 

0.02 

-0.03 

-0.02 

0.58 

1.71 

1.93 

1.98 

2.00 

£o/2" 

7.11E-3 

3.30E-3 

4.07E-3 

4.43E-3 

4.53E-3 

4.56E-3 

3.05E-3 

9.32E-4 

2.45E-4 

order 

- 

1.11 

-0.30 

-0.12 

-0.03 

-0.01 

0.58 

1.71 

1.93 

eo/2^ 

7.19E-3 

1.99E-3 

5.10E-4 

6.84E-4 

1.02E-3 

l.lOE-3 

1.13E-3 

1.14E-3 

7.61E-4 

order 

- 

1.85 

1.96 

-0.42 

-0.58 

-0.11 

-0.04 

-0.01 

0.58 

eo/2'' 

7.07E-3 

1.70E-3 

4.49E-4 

2.61E-4 

8.81E-5 

1.68E-4 

2.54E-4 

2.77E-4 

2.83E-4 

order 

- 

2.06 

1.92 

0.78 

1.57 

-0.93 

-0.60 

-0.13 

-0.03 


7.05E-3 

1.71E-3 

4.23E-4 

1.09E-4 

3.91E-5 

6.01E-5 

2.18E-5 

4.20E-5 

6.35E-5 

order 

- 

2.04 

2.02 

1.96 

1.48 

-0.62 

1.46 

-0.95 

-0.60 

eo/2^ 

7.05E-3 

1.71E-3 

4.22E-4 

1.05E-4 

2.61E-5 

1.37E-5 

6.98E-6 

1.50E-5 

5.48E-6 

order 

- 

2.04 

2.03 

2.01 

2.01 

0.93 

0.97 

-1.10 

1.45 


To obtain a second order limiting Pauli-type equation model, we formally drop the small compo¬ 
nents in (I2.10I) - (I2.11I) to get 

(4.5) a;) = x) -I- 4'p(t, x) -I- O(e^), 

where 4'e := 4'e(t,x) G and 4'p := 4'p(t,x) G satisfy the Pauli-type equations 

(4.6) + n+ (tp^e) , -f H. (IP^p) , X G M, t > 0, 

with V = \/I — e^A — / and initial value as 


(4.7) 'he(0, x) = n+<i)(0, x), 4'p(0, x) = n_$(0, x), X G M. 

To investigate numerically convergence rates of the above limiting models (14.31) and (j4.5l) to the Dirac 
equation, we solve numerically the Schrodinger equation (14.41) to obtain and the Pauli-type 

equation (14.61) to get (4'e,4'p), by the TSFP method and the EWI-FP method m, respectively. 
The solution $ of the Dirac equation (11.81) is computed by the MTI-FP method and we can study 
convergence rates of Dirac equation (11.81) to (14.31) and (14.51) . respectively. All the computations are 
done on the bounded interval D = (—128,128) with fine mesh h = 1/16 and time step r = 10“"^, 
which are sufficiently small such that the numerical error can be neglected. In order to quantify the 
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Fig. 4.1. Error functions F'sch(^) Epauit) for different e. 


convergence, we introdnce the error functions 


Eschit) = •)-£ •)e2 


L2 


£’pau(i) = 


L2 


i > 0. 


Fig. 14.II depicts the evolution of the errors Esch{t) and Epau{t), and we can conclude that the limiting 
model of the Schrodinger equation (14.31) is linearly accurate at 0{e), while the limiting model of the 
Pauli-type equation (14.51) is quadratically accurate at O(e^). In particular, both the errors Esch{t) and 
Epanit) are observed to grow linearly in time, i.e. 

(4.8) Eschit) < (Cl + C2t)e, Epauit) < {C 3 + C4^)e^ t > 0, 


where Ci, C 2 , C 3 and are positive constants independent of time t > 0 and e G (0,1]. We find that 
(14.61) is the same second order approximate limiting model as the Pauli equation [3TJ[33] for the Dirac 
equation (11.81) in the nonrelativistic limit regime. 
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5. Conclusion. A multiscale time integrator Fourier pseudospectral (MTI-FP) method was pro¬ 
posed and rigorously analyzed for the Dirac equation involving a dimensionless parameter e € (0,1], 
which is inversely proportional to the speed of light. The main difficulty of the problem is that the 
solution highly oscillates with O(e^) wavelength in time when 0 < e <C 1. The key ideas in designing 
the MTI-FP method included a proper multiscale decomposition of the Dirac equation and the use of 
the Gautschi type exponential wave integrator in time discretization. Rigorous error analysis showed 
that the MTI-FP method is uniformly convergent in spatial discretization with spectral accuracy, and 
uniformly convergent in temporal discretization with linear order for e G (0,1], while the temporal 
accuracy is optimal with quadratic convergence rate when either e = 0(1) ore < t. This result 
significantly improves the error bounds of the existing numerical methods for the Dirac equation in 
the nonrelativistic limit regime. Numerical results confirmed the error estimates and suggested our 
error bounds are sharp and optimal. Convergence rates of the Dirac equation to its limiting first order 
Schrodinger equation model and second order Pauli-type equation model were observed numerically. 
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